#############################
##CPS: Politics of Pessimism
##Replication file: Maps
##AX
##############################

########################################################
## Setting working directory
########################################################

library(foreign)
library(descr)
library(plotrix)
library(readstata13)
library(sp)
library(plyr)
library(sf)
library(raster)
library(descr)
library(RColorBrewer)
library(plotrix)
library(ggplot2)
library(classInt)
library(sp)
library(tmap)
library(tidyverse)

##Read in aspirations data
setwd("~/Dropbox/Trump Polls/CPS Replication/Replication_final/Data/")
aspirations <- read.dta13(file = "gallup_prelim_clean2015_collapsedbyfipscode.dta")

##Read in U.S. TIGER state shapefiles
setwd("~/Dropbox/Trump Polls/CPS Replication/Replication_final/Data/tl_2016_us_county/")
US.states  <- shapefile('tl_2016_us_county.shp')
head(US.states@data)

# returns string w/o leading or trailing whitespace
trim <- function (x) gsub("^\\s+|\\s+$", "", x)

US.states@data$STATEFP <- trim(US.states@data$STATEFP)
US.states@data$COUNTYFP <- trim(US.states@data$COUNTYFP)

US.states@data$fipscode <- paste(US.states@data$STATEFP,US.states@data$COUNTYFP, sep ='')
head(US.states@data)

########################################################
## Trim sample to remove Guam, outer areas
########################################################

US.states <- subset(US.states, STATEFP != "03")
US.states <- subset(US.states, STATEFP != "07")
US.states <- subset(US.states, STATEFP != "14")
US.states <- subset(US.states, STATEFP != "43")
US.states <- subset(US.states, STATEFP != "72")
US.states <- subset(US.states, STATEFP != "52")
US.states <- subset(US.states, STATEFP != "02")
US.states <- subset(US.states, STATEFP != "15")

########################################################
## Merge and process data
########################################################

# Merge data
US.states@data$order <- 1:nrow(US.states@data)
US.states@data <- merge(US.states@data, aspirations, by = "fipscode", all.x=T)
head(US.states@data)
summary(US.states@data$neg_ladderfuture)
# ================================================================
# plot
# ================================================================
#brewer.pal.info
myPaletteRed <- brewer.pal(9,'Reds')
createColors <- colorRampPalette(myPaletteRed)
myColorsRed <- createColors(nrow(US.states@data))

## use quantile to create cut points in dependent variable
US.states@data <- US.states@data[order(US.states@data$neg_ladderfuture),]
numberColors <- length(unique(myColorsRed))

cuts <- quantile(US.states@data$neg_ladderfuture, probs = (0:numberColors)/numberColors, na.rm = T) 
US.states@data$Ecat <- cut(US.states@data$neg_ladderfuture, breaks = unique(cuts), include.lowest=T)

# merging Ecat and its colors
col2 <-"#7f1701" # Color for all extreme right colors
colorOfCats <- data.frame(Qcat=levels(US.states@data$Ecat), colors=c(unique(myColorsRed)[1:395], col2))
US.states@data <- merge(US.states@data,colorOfCats, by.x='Ecat', by.y="Qcat", all.x=T)

# going back to the original order of the data in the shapefile
US.states@data <- US.states@data[order(US.states@data$order), ]

pdf("covidtrump_mapaspirations.pdf", width=7, height=8)
plot(US.states, col=as.character(US.states@data$colors), 
     lty=0, main= 'Reduced Aspirations, 2015')
plot(US.states, add=T) 

# legend
color.legend(xl=-72, xr=-64,yb=-26,yt=-27, 
             legend=c('0%','','100%'), 
             rect.col=myColorsRed, gradient='x',
             cex=.7, pos=c(1,1,1))
#color.legend(xl=-66, xr=-64, yb=-28, yt=-29,
#             legend=('> 61.1%'), rect.col=col2, cex=.7, pos=c(1,1,1))
text(x=-68,y=-24,label='Future Aspirations, (Pct. Point)', cex=.8)
dev.off()


